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A method for manipulating phonocardiography signals. 

5 The invention relates to a method for assisting medically trained persons in 

diagnosing and classifying medical conditions of the heart by means of auscultation. 

Signals obtained by means of a transducer are phonocardiography representations of 
sounds traditionally listened to by means of a stethoscope. Training in auscultation 

1 0 takes a long time and requires an aptitude for recognising and classifying aural cues, 
frequently in a noisy environment 20-30 different conditions may need to be 
differentiated, and within each, the severity evaluated. Furthermore, there may be 
combinations among these. These factors contribute to explaining why not all 
physicians perform equally well when diagnosing heart conditions, and why it may 

1 5 be time-consuming. 

Use of modern signal processing has been proposed to enhance the relevant aural 
cues, and in some cases it has been demonstrated that a transposition into a different 
frequency range or repetition of a repetitive signal at a different repetition frequency 
20 from that originally recorded may be advantageous. Among the patents that describe 
this kind of signal processing the following may be mentioned: 

US 4,220,160 uses the heart sound picked up by a microphone to modulate a higher 
and more audible carrier frequency. US 4,424,815 is a very fundamental patent 

25 employing a range of solutions, from magnetic drum (i.e. a non-digital procedure) to 
a tapped delay line for stretching time while conserving frequency content. US 
4,528,689 describes a time stretching digital stethoscope that records several cycles, 
and obtains a "slow-down factor" proportional to the number of cycles recorded. 
Each part of a sound between zero crossings is repeated with opposite polarity, 

30 thereby obtaining a signal that sounds like original signals but lasting longer. Only 
whole-number factors are possible. US 4,792,145 describes how frequency scaling is 
obtained by multiplying each Fourier spectral component by a factor. The original 
relationship between frequency components is retained, and it is purported that the 
audibility is improved. Similarly, the time may be scaled. 
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All of the above solutions have severe limitations, either because of an increase in 
non-linear distortion and noise, or because the factor for stretching can only be 
whole-numbered. 

5 

According to the invention, these disadvantages may be avoided by synthesizing a 
"clean" signal consisting of a sum of sinusoids and to perform a transformation on 
these sinusoids in order that the relative amplitudes are maintained while the 
frequency range or time axis are changed in a continuous range of scales. In a further 
1 0 refinement of this method, the scale may be linked to other phenomena by means of 
an auto-scale function. 

The invention will be more fully described in the following with reference to the 
drawing, in which 

15 

Fig. 1 shows a block diagram of the operations performed on an input signal in order 
that a sinusoidal model is obtained for the signal, which is added to obtain an 
output signal. 

Fig. 2 shows the procedure of block segmentation, 
20 Fig. 3 shows the estimation of a finite number of spectral peaks, 

Fig. 4 shows results of the peaks search 

Fig. 5 shows the matching procedure 

Fig. 6 shows the final results of the peak matching 

Fig. 7 shows boundary condition for the interpolation functions, 
25 Fig. 8 shows interpolation of instantaneous phase function and amplitude function 

Fig. 9 shows a spectrogram of the original heart sound, 

Fig. 10 shows a spectrogram of the sinusoidal coded heart sound, 

Fig. 1 1 shows a Time-Frequency scaling procedure 

Fig. 12 shows a spectrogram of a frequency shifted heart sound (p=l,^=l,r=100Hz), 
30 Fig. 13 shows a spectrogram of a frequency scaled heart sound (p=l ,^=2,r=0), 
Fig. 14 shows a spectrogram of a non-scaled sinusoidal coded heart sound 
07=l,^=l,r=0), 

Fig. 1 5 shows a spectrogram of a joint time-frequency scaled heart sound (p as 2,^= s 2, 
a-O), 
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Fig. 16 shows a spectrogram of a time scaled heart sound (p=2,q=l,r=0), and 
Fig. 17 shows a spectrogram of the original heart sound (p=l,q=l,r=0). 

The mathematical foundation for this procedure is described in the following 
5 paragraphs. 

The sinusoidal model was originally developed for speech applications by McAulay 
and Quantieri and a thorough description can be found in Robert J. McAulay and 
Thomas F. Quatieri, "Speech Analysis/Synthesis Based On A Sinusoidal Representa- 

1 0 tion", IEEE Transactions on Acoustics, Speech, and Signal Processing, Vol. ASSP- 
34, NO. 4, August 1986, page 744-754, and Thomas F. Quatieri and Robert J. 
McAulay, "Speech Transformations Based On A Sinusoidal Representation", IEEE 
Transactions on Acoustics, Speech, and Signal Processing, Vol. ASSP-34, NO. 6, 
December 1986, page 1449-1464. Because the spectral structure of a speech signal 

1 5 and a heart signal are very different (a heart beat does not have pitch like voiced 
speech, nor a broadband structure like unvoiced speech), the original sinusoidal 
model will not be able to model a heart beat properly. A sinusoidal model specific 
for heart sounds is developed by exploring the spectral characteristics of the heart 
beat, when tracking spectral components in the signal. 

20 

The sinusoidal model is based upon sinus functions. The synthesis part of the model 
consists in adding a group of sinusoidal functions: 

K0 - Z ^<<>cos(fi A (/)) 

25 The number of functions L(t) will change as a function of time, depending on the 

complexity of the signal. Each sinusoidal function (frequency track) is controlled by 
the total phase function Q(/), and the amplitude L(t). The main objectives of the 
analysis part of the model is to establish the content of these two functions. 

A fundamental relation for the sinusoidal model is the relationship between total 
30 phase and the instantaneous frequency: 



Eq. (2): 

(o(/) - JjQ(0 [rad/s] 



Substituting the expression for the total phase into eq. (1), the final expression for the 
sinusoidal model is obtained: 



5 Eq. (3): 



10 




The structure of a traditional sinusoidal model is illustrated in Fig. 1. The input 
signal *(/) is obtained from a sensor for heart sounds placed in contact with the body 
under examination. The output signal y(f) is made available to the medically trained 
person via e.g. headphones after suitable amplification. Furthermore it may be used 
for a visual display and comparison of frequency spectra on a display screen. 



IS 



The model may be divided into three parts: an Analysis part where the signal 
parameters are measured, an Interpolation and transformation part, where Q(/) and 
A(t) are interpolated from one window position to the next window position, and a 
Synthesis part where the transformed heart signal is reconstructed. 

The content of each block of Fig. 1 will be elaborated more in the following sections. 



The first step of the Analysis is to measure the frequency content of the signal. As 
the frequency content is changing as a function of time, a Short Time Fourier 
20 Transform (STFT) is used to track the frequency content. The STFT segments the 
signal into overlapping segments (75%), and each segment is Fourier transformed 
using a Discrete Fourier Transform (DFT). The segmentation procedure is illustrated 
in Fig. 2. 



25 When a signal segment has been Fourier transformed, the most dominant frequency 
component in the spectrum has to be found, and used as local "check points" for a 
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frequency track over time. The original sinusoidal model finds every peak in the 
spectrum, and associates a sinusoidal generator to each peak. In the case of a voiced 
speech signal, the signal is periodic and the corresponding spectrum will be discrete 
(pitch). As the spectrum is discrete, the known peak search algorithm will find a 
5 large number of peaks. However, a heart beat is non-periodic and has a nature like a 
Gaussian function modulated by a low-frequency sinusoid. The corresponding 
spectrum is also Gaussian, and the known peak search algorithm will model this 
heart beat by only one sinusoidal function, where the amplitude is equal to the peak 
level. This will not model the heart beat correctly, compared to the remaining 
1 0 spectrum. However, if the amplitude of the sinusoidal function instead is set equal to 
the area below the curve defining the "hill" corresponding to the peak, the heart beat 
will be properly modelled. The peak searching procedure is illustrated in Fig. 3. 

First the index for the maximum peak is found. Secondly, the area of the "hill" is 
1 5 estimated, and the complete "hill" is deleted from the spectrum. The procedure is 
repeated with the modified spectrum, until a maximum number of peaks has been 
found, or the peak level is below a pre-defined threshold. For certain heart sounds as 
many as 50 peaks corresponding to 50 sinusoids may be found. 

20 The result from the peak search is a table of peaks for each position of the window 
function. The table corresponding to window position m will be termed spec m. The 
distance between each window position is T seconds. This is illustrated in Fig. 4. 

The peaks in spec m now have to be matched together with peaks in spec m+1 , in 
25 order to construct a frequency track that the total phase function must follow. 
Many different solutions for this matching problem have been proposed in the 
literature, as no optimal solution exits. The matching strategy used in the present 
invention is based on the following three rules: 

- high energy peaks are matched first 
30 - tracks are not allowed to cross. 

- the frequency differences between two matched peaks must be below certain 
limits. 
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First, the algorithm finds the maximum peak in spec m. This peak has to be matched 
with a peak in spec m+L For this purpose, the algorithm searches spec m+1 for the 
maximum peak within lower and upper limits. These limits are controlled by two 
factors. The distance to the nearest peak in spec m and the maximum allowed 
5 frequency change during T seconds. Before a match is accepted, it must be verified, 
i.e. it must not cross a previously found match. When the final match is verified, the 
two peaks are no longer visible to the matching algorithm. Unmatched peaks in spec 
m are declared dead, and unmatched peaks in spec m+1 are declared born. The 
procedure is illustrated in Fig. 5, and an example of the matching process is 
10 illustrated in Fig. 6. 

The synthesis will be designed for T seconds at a time, that is from window position 
m to window position m+1, and this T second reconstructed signal will be termed 
frame m. The number of born, continued, and dead tracks in frame m will be termed 
15 L\(m) 9 L2(m), and L3(m). For each frame m, the born, continued and dead tracks are 
synthesized separately and added to form the final signal: 

Eq. (4): 

AO -*(0+*(0+j£(0 

L\{m) L2(m) £3(m) 

= £ (0003(^(0)+ £ 4a(0cai(n£(0)+ £ /^(0 cos (0^(0) 
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Each track k in frame m has to be synthesised by a sinusoidal function 
i4*(0cos(Q*(0« The track is up to now only specified by the start and end point. This 
is illustrated in Fig. 7. 

25 A m k is the amplitude of the Fourier coefficients, co m * is the frequency of the Fourier 
coefficient and y m k is the phase offset for the current track. The phase offset is 
increased by an average phase step dependent on the start and stop frequency of the 
track (see Fig. 8). The original model is based on the phase of the Fourier 
coefficients. The amplitude is interpolated by A k (t\ and the instantaneous phase is 

30 interpolated by Q/c(0- The interpolation functions are only constrained by the start 
and stop conditions: 
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Eq. (5): 





- A m 
~ A k 


A k {T) = AS* 1 

fn m + l 

m CO. + CD. 


^(0) 




*V(0) 







The amplitude function is not critical, so a linear interpolation function will be used 
for /i*(0- However, the instantaneous phase function is more critical because 
5 any discontinuity in the phase progress will deteriorate the perceptual quality of the 
final sound. In the following, the synthesis of born, continued, and dead tracks will 
be designed separately. 

In the case of a continued track, a second order polynomial will be used to 
1 0 interpolate a smooth phase progress during a frame: 

Eq.(6): 

Q(/) = a l +a 2 / + a 3 -/ 2 
This is illustrated in Fig. 8. 

15 

The coefficients are determined by inserting the conditions from Eq, (5) into Eq. (6): 
Eq. (7): 



fi(0) = 




- ,„ m 

- <P* 


Q'(0) = 


a i 


m 

= <o A 


Q(7} = 


a l 


m m + l m + 1 m 

+a 2 -r+o 3 -r - cp"+r- 2 * — =>a 3 - 2T 



The linear interpolation of the amplitude function is achieved by: 
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Eq. (8): 



The final equations for the synthesis of a continued track k in frame m is given by: 



Continued track 
5 Eq. (9): 



fl* + 1 m 

Q *(0 = <P* +G H ' + J? 



Two special cases remain to be considered. When a track is declared dead or 6o/yi, 
we do not have both start and end conditions for them. 

1 0 In the case of a born track, an instantaneous phase function will be initiated from the 
begining of frame k, with A m u = 0, oA « ©V <p m * = 0, and the second order 
instantaneous phase function is reduced to a linear function: 
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Born track: 
Eq. (10): 



n*<0 




4Te> 


,w + 1 



When a track is declared dead, the track is terminated at the end of frame m using 
A m k ' = 0, a/V ' = <»A, and q> V 1 =T- <p m *, and the second order instantaneous 
phase function is again reduced to a linear function: 



20 



Dead track: 
Eq. (11): 



<(0 = 



In order to explain more fully the sinusoidal model according to the invention it is 
prepared for implementation in a mathematical simulation program, such as MatLab 
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(TM), and the interpolation functions are hence converted to discrete index. The 
discrete time version is obtained using the following substitutions: 

Eq. (12): 



/ = 


n • A 




T = 


N- A 




m 


m 




<B k = 


N DFT 


A 



These substitutions are now inserted into the main equations for the sinusoidal 



The time index n is running from 0 to N. The discrete amplitude function is given by: 
Eq.(l4): 



model. 



Eq.(13): 




A 

A^[n] = 4?(»-A) = ^-/i-A--*- 




The final equations for a discrete-time continued track is: 



Discrete continued track: 
Eq.(l5): 



nn»] = <pr+«vr-r22-+i. i -(vr 1 -v?) 



N-N, 



DFT 




10 

The discrete function for a born and dead track can likewise be found following the 
same approach: 



Discrete born track: 
Eq. (16): 



m + 1 



n k [n] = n • A • — — « * • v k • — 



*DFT 



Discrete dead track: 
Eq. (17): 



10 The result of coding a heart sound with severe systolic murmurs is illustrated in Figs. 
9 and 10. 



Using the novel sinusoidal model developed above, we are now in a position to make 
modifications to the original signal, which would normally be very difficult to 
1 5 obtain. As the model contains a very precise description of how the frequency 

content of the signal is changing over time, it provides a good foundation to make 
signal transformations like time stretching (without changing the frequency 
structure), frequency stretching (without changing the time progress), and joint time 
frequency stretching. 

20 

The sinusoidal model may be extended to a time-frequency scaled sinusoidal model 
with small modifications. The joint time frequency scaling is obtained by mapping 
track information at time t and frequency co in the original model to a new time 
position r « p • /, where p is the time scaling factor, and new frequency position aT = 
25 q • /, where q is the frequency scaling factor. A pure frequency scaled model, may be 
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obtained by setting p = 1 , and a pure time scaling model may be obtain by setting q = 
1. The procedure is illustrated in Fig. 1 1. 

The time scaling is implemented in the sinusoidal model by changing the frame T to 
p • T in the interpolation functions and A(t). Frequency scaling can be 
5 implemented by shifting the peaks found by the peak searching algorithm from a> to 
q • cd: 



T = 


P 


■T 


» = 


q 





Eq.(18): 

When q • co crosses half the sampling frequency, the track is removed. Every time a 
I o born or continued track has been synthesized, the phase offset (pV ; has to be 
updated to: 

Eq.(19): 



Other frequency modifications will also be possible. Heterodyne modification is 
1 5 obtained by shifting peaks from to to co + r, where r is the amount of the spectrum 
movement. 

The time-scaled sinusoidal model may form the foundation for an autoscale 
algorithm, which scales arbitrary heart sounds to 60 beats per minute (bpm). This 
20 autoscale algorithm is obtained using a scale factor p: 
Eq. (20): 

where h is the heart rate. There is no problem in following a changing heart rate 
25 dyamically. 



# 
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In Figs. 12 to 17, the results of applying the above described procedure are shown 
when used on a heart sound. 



It will be understood that once the signal has been converted to digital representation 
5 of data, its manipulation may take place in dedicated processors, RISC processors or 
general purpose computers, the outcome of the manipulation being solely dependent 
on the instructions performed on the data under the control of the program written for 
the processor in order to obtain the function. The physical location of the data at any 
one instant (i.e. in varying degrees of processing) may or may not be related to a 
1 0 particular block in the block diagram, but the representation of the invention in the 
form of interconnected functional blocks provides the skilled person with sufficient 
information to obtain the advantages of the invention. 

The foregoing description of the specific embodiments will so fully reveal the 
1 5 general nature of the present invention that others skilled in the art can, by applying 
current knowledge, readily modify or adapt for various applications such specific 
embodiments without undue experimentation and without departing from the generic 
concept, and therefore, such adaptations and modifications should and are intended 
to be comprehended within the meaning and range of equivalents of the disclosed 
20 embodiments. It is to be understood that the phraseology or terminology employed 
herein is for the purpose of description and not of limitation. The means, materials, 
and steps for carrying out various disclosed functions may take a variety of forms 
without departing from the invention. 

25 Thus, the expressions "means to ... " and "means for or any method step 

language, as may be found in the specification above and/or in the claims below, 
followed by a functional statement, are intended to define and cover whatever 
structural, physical, chemical, or electrical element or structure, or whatever method 
step, which may now or in the future exist which carries out the recited functions, 

30 whether or not precisely equivalent to the embodiment or embodiments disclosed in 
the specification above, i.e., other means or steps for carrying out the same function 
can be used; and it is intended that such expressions be given their broadest 
interpretation. 
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PATENT CLAIMS 

1 . A method for manipulating phonocardiography signals obtained from a 
transducer and subjected to digital spectral analysis, characterised in that 
5 the signals are converted on a running basis to a sinusoidal model representation, that 
a time axis scaling is defined and used to control the amplitudes and phases of the 
sinusoids, which are subsequently added to create a time scaled representation of the 
phonocardiographic signal. 

10 2. A method for manipulating phonocardiographic signals obtained from a 

transducer and subjected to digital spectral analysis, characterised in that 
the signals are converted on a running basis to a sinusoidal model representation, that 
a frequency axis scaling is defined and used to control the amplitudes and phases of 
the sinusoids, which are subsequently added to create a frequency scaled 

1 5 representation of the phonocardiographic signal. 

3. A method according to claim 1, characterised in comprising the 
following steps: 

- obtaining the frequency content of the signal by applying a Short Time Fourier 
20 Transform with overlapping segments, a Discrete Fourier Transform being 

performed on each segment, 

- performing a frequency peak search on each segment by repeated removal of the 
highest spectral components, identifying each peak by its frequency value 

- repeating the peak search until a maximum number of peaks have been identified 
25 or until the peak level of the last peak is below a preset minimum, 

- establishing a segment-by-segment map of spectral peaks, said peaks forming a 
track over time 

- adjusting the phase of sine generators centered on the frequencies of the tracks 

- summing the outputs of all sine generators active at any one instant. 

30 

4. A method according to claim 2, characterised in comprising the 
following steps: 

- obtaining the frequency content of the signal by applying a Short Time Fourier 
Transform with overlapping segments, a Discrete Fourier Transform being 
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performed on each segment, 

performing a frequency peak search on each segment by repeated removal of the 
highest spectral "hills" identifying each peak by its frequency value 
repeating the peak search until a maximum number of peaks have been identified 
or until the peak level of the last peak is below a preset minimum, 
establishing a segment-by-segment map of spectral peaks, said peaks forming a 
track over time 

adjusting the phase of sine generators centered on the frequencies of the tracks 
summing the outputs of all sine generators active at any one instant. 
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ABSTRACT 



A heart signal comprising both first and second heart sounds and murmurs occurring 
5 due to various heart conditions is modelled by a sum of sinusoids automatically 
selected to represent the heart sound without a noise component. In an analysis step 
the signal parameters are measured, in an interpolation and transformation part Q(f) 
and A(t) are interpolated from one window position to the next window position, and 
a synthesis part reconstructs the transformed heart signal. The transformations are 
10 transformations in either time scale (with unchanged frequency content) or frequency 
scale (with unchanged time progression), or both. Furthermore a dynamic 
transformation of an arbitrary heart signal to a fixed number of beats per minute is 
made possible. 
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